require(Rcpp)
options(rasterTmpDir='C:/temp/')
gnut_fn <- file.path(gnut_dir,'spam2010V2r0_global_P_GROU_A.tif')

output_version <- '2020-10-20'

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01') %>%
  st_set_crs(.,4326)

ntl.zones.prj <- ntl.zones %>%
  # 102022
  st_transform(.,crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

# convert polygons to points 
ntl.zones.pts = st_as_sf(as.data.frame(M.ntl.zones), coords = c("X", "Y"), crs = 4326)
ntl.zones = st_drop_geometry(subset(ntl.zones,select=c(OBJECTID, GID_0)))
ntl.zones.pts = dplyr::bind_cols(ntl.zones.pts,ntl.zones)

##
## GROUNDNUT : SPAM 2010v2
##
ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01')

  gnut.global <- raster(gnut_fn)
  
  # manually set outside for 5x5
  ntl.extent <- c(-5,30,-5,30)
  gnut <- crop(gnut.global,ntl.extent,snap='out')

## - create measure avg same cell, avg 5 cells from cell, avg 10 cells from cell
gnut.diam.crop <-  sqrt(raster::area(gnut))
decdeg050km = 50 / cellStats(gnut.diam.crop,'mean')

lendd050km = res(gnut.diam.crop)[1] * decdeg050km  # afclima.crop??

decdeg100km = 100 / cellStats(gnut.diam.crop,'mean')
# 1deg = 111 km  (or 60 nautical miles)
lendd100km = res(gnut.diam.crop)[1] * decdeg100km

fw5x5 <- focalWeight(gnut, lendd050km, type='circle')
fw10x10 <- focalWeight(gnut, lendd100km, type='circle')

gnut050 <- focal(gnut, w=fw5x5, fun=sum,na.rm=TRUE)
gnut100 <- focal(gnut, w=fw10x10, fun=sum,na.rm=TRUE)


# sum production in cell
ntl.zones$gnut <- exactextractr::exact_extract(gnut,ntl.zones,fun='sum')
ntl.zones$gnut050 <- exactextractr::exact_extract(gnut050,ntl.zones,fun='sum')
ntl.zones$gnut100 <- exactextractr::exact_extract(gnut100,ntl.zones,fun='sum')

save.dta13(as.data.frame(st_drop_geometry(subset(ntl.zones,select=c(OBJECTID,GID_0,gnut,gnut050,gnut100)))), 
           paste0(wkdir,"/proc_data/GNUT",output_version,".dta"))
